additional_analyses
Data preparation
Data summary
Traits correlation
- remove the traits that are too correlated?
Code
M <-cor(numeric_traits[, c(-1)])
ggcorrplot::ggcorrplot(M, hc.order = TRUE, type = "lower",
lab = TRUE, tl.cex = 9, lab_size = 3)
- Traits types
CWM bootstrap
By depth layer
each trawl is a replica
Utilisation traitstrap R package
Bootstrap non parametric, échantillonnage dans les proportions de la biomasse des espèces
pour le remplissage des valeurs manquantes : si manque données pour 1 individus alors va chercher la valeurs sur un ind de la même profondeur, si pas possible alors ind à une autre profondeur mais de la même couche et si impossible ind d’une autre couche de profondeur
Code
spxcom.matrix <- depth_fish_biomass %>%
t() %>%
as.data.frame() %>%
relocate("Epipelagic", "Upper mesopelagic", "Lower mesopelagic","Bathypelagic" ) %>%
tibble::rownames_to_column("species") %>%
arrange(species) %>%
tibble::column_to_rownames("species") %>%
as.matrix()
spxtraits.matrix <- fish_traits %>%
mutate(across(16:23, ~ case_when(. == "P" ~ 1,
. == "A" ~ 0,
TRUE ~ as.numeric(.))),
across(24, ~ case_when(. == "A" ~ 1,
. == "B" ~ 2,
. == "C" ~ 3,
TRUE ~ as.numeric(.)))) %>%
select(-c(gill_raker_types, oral_gape_axis)) %>%
as.matrix()
# Remove the "[,1]" suffix from column names
names(spxtraits.matrix) <- gsub("[,1]", "", names(spxtraits.matrix))
#check rownames
#rownames(spxtraits.matrix) == rownames(spxcom.matrix)
result_CWM <- FD::functcomp(spxtraits.matrix, t(spxcom.matrix))
#FD::functcomp(spxtraits.matrix, t(spxcom.matrix), CWM.type = "all")
# Calculate Total biomass
total_biomass <- colSums(spxcom.matrix)
# Calculate Relative biomass
sp_rel_biomass <- t(spxcom.matrix) / total_biomass
# Transpose the Relative biomass Matrix for Display
t_sp_rel_biomass <- t(sp_rel_biomass)
total_sum <- colSums(t(sp_rel_biomass))
# Initialize an empty data frame to store results
CWM_df <- data.frame(
depth_layer = character(),
trait = character(),
total_sum = numeric(),
weighted_mean = numeric(),
stringsAsFactors = FALSE
)
library(traitstrap)
# Trait
trait_boot <- morpho_data%>%
inner_join(metadata) %>%
select(-c(individual_code, years, longitude_start,
latitude_start, longitude_end, longitude_end,
volume_filtered, distance)) %>%
mutate(
eye_size = eye_diameter / head_depth,
orbital_length = eye_diameter / standard_length,
oral_gape_surface = mouth_width * mouth_depth / body_width * body_depth,
oral_gape_shape = mouth_depth / mouth_width,
oral_gape_position = distance_upper_jaws_bottom_head / head_depth,
lower_jaw_length = lower_jaw_length / standard_length,
head_length = head_length / standard_length,
body_depth = body_depth / standard_length,
pectoral_fin_position = distance_pectoral_bottom_body / body_depth_pectoral_insertion,
pectoral_fin_insertion = prepectoral_length / standard_length,
transversal_shape = body_depth / body_width,
dorsal_fin_insertion = predorsal_length / standard_length,
eye_position = eye_height / head_depth,
operculum_volume = operculum_depth / operculum_width,
gill_outflow = operculum_width,
caudal_throttle_width = caudal_peduncle_min_depth
) %>%
select(
depth,
species,
eye_size,
orbital_length,
gill_outflow,
oral_gape_surface,
oral_gape_shape,
oral_gape_position,
lower_jaw_length,
head_length,
body_depth,
pectoral_fin_position,
pectoral_fin_insertion,
transversal_shape,
caudal_throttle_width,
dorsal_fin_insertion,
eye_position,
operculum_volume,
ventral_photophores,
gland_head,
chin_barbel,
small_teeth,
large_teeth,
fang_teeth,
retractable_teeth,
internal_teeth
) %>%
mutate_at(vars(ventral_photophores,
gland_head,
chin_barbel,
small_teeth,
large_teeth,
fang_teeth,
retractable_teeth,
internal_teeth),
funs(ifelse(. == "P", 1, ifelse(. == "A", 0, .)))) %>%
mutate(across(all_of(c("ventral_photophores",
"gland_head",
"chin_barbel",
"small_teeth",
"large_teeth",
"fang_teeth",
"retractable_teeth",
"internal_teeth")), as.numeric)) %>%
tidyr::pivot_longer(!c(species,depth), names_to = "trait", values_to = "values")%>%
# add column with depth layer
mutate(
depth_layer = case_when(
between(depth, 0, 174) ~ "Epipelagic",
between(depth, 175, 699) ~ "Upper mesopelagic",
between(depth, 700, 999) ~ "Lower mesopelagic",
between(depth, 1000, 2000) ~ "Bathypelagic"))
# Community
community <- rbind(data_biomass_2002_2019, data_biomass_2021, data_biomass_2022)%>%
as.data.frame()%>%
rename("species"="Nom_Scientifique",
"station"="Code_Station") %>%
left_join(metadata) %>%
select(species, Tot_V_HV, depth, volume_filtered)%>%
# add column with depth layer
mutate(
depth_layer = case_when(
between(depth, 0, 174) ~ "Epipelagic",
between(depth, 175, 699) ~ "Upper mesopelagic",
between(depth, 700, 999) ~ "Lower mesopelagic",
between(depth, 1000, 2000) ~ "Bathypelagic"))%>%
replace(is.na(.), 0)%>%
group_by(species, depth)%>%
mutate(biomass=sum(Tot_V_HV))%>%
select(-c(Tot_V_HV))%>%
distinct()%>%
select(-c(volume_filtered)) %>%
filter(biomass>0) %>%
mutate(species = case_when(
species == "Nannobrachium_atrum"~"Lampanyctus_ater",
species == "Cyclothone"~"Cyclothone_sp",
species == "Stomias_boa_boa"~"Stomias_boa",
TRUE ~ species
))
trait_filling <- trait_fill(
# input data (mandatory)
comm = community,
traits = trait_boot,
# specifies columns in your data (mandatory)
abundance_col = "biomass",
taxon_col = "species",
trait_col = "trait",
value_col = "values",
# specifies sampling hierarchy
scale_hierarchy = c("depth_layer", "depth"),
# min number of samples
min_n_in_sample = 5
)
# run nonparametric bootstrapping
np_bootstrapped_moments <- trait_np_bootstrap(
trait_filling,
nrep = 100
)
np_bootstrapped_moments_plot <- np_bootstrapped_moments %>%
mutate(trait= gsub("_"," ", trait))%>%
filter(
trait %in% c(
"caudal throttle width",
"oral gape surface",
"gill outflow",
"large teeth",
"small teeth",
"orbital length",
"eye size",
"transversal shape",
"operculum volume",
"ventral photophores",
"internal teeth",
"eye position",
"oral gape shape"
)
)
np_bootstrapped_moments_plot$trait <- factor(
np_bootstrapped_moments_plot$trait,
levels = c(
"caudal throttle width",
"oral gape surface",
"gill outflow",
"large teeth",
"small teeth",
"orbital length",
"eye size",
"transversal shape",
"operculum volume",
"ventral photophores",
"internal teeth",
"eye position",
"oral gape shape"
)
)
np_bootstrapped_moments_plot$depth_layer <- factor(np_bootstrapped_moments_plot$depth_layer,
levels = c("Epipelagic", "Upper mesopelagic", "Lower mesopelagic", "Bathypelagic"))
# Mean
ggplot(np_bootstrapped_moments_plot, aes(x=depth_layer, y=mean)) +
geom_boxplot(aes(col=depth_layer, fill=depth_layer), alpha=0.08) +
scale_fill_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A"))+
scale_color_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A"))+
guides(col="none", fill="none")+
labs(y="Bootstrapped CWM")+
facet_wrap(~trait, scales = "free")+
theme_light()+
theme(axis.text.x = element_blank(),
axis.title.x = element_blank(),
panel.grid.minor = element_blank(),
axis.title.y.left = element_text(size =14),
strip.text = element_text( face="bold"),
axis.title.y = element_text(size=15),
axis.text.y = element_text(size=14),
strip.background=element_rect(fill="white"),
strip.text.x = element_text(size = 14, color = "black"))
Code
ggsave("CWM_boot.png", path = "figures/additional_analyses", dpi = 700, height = 10, width = 12)PCA CWM
Code
sum_boot_moment <- trait_summarise_boot_moments(
np_bootstrapped_moments
) %>%
ungroup() %>%
filter(
trait %in% c(
"caudal_throttle width",
"oral_gape_surface",
"gill_outflow",
"large_teeth",
"small_teeth",
"orbital_length",
"eye_size",
"transversal_shape",
"operculum_volume",
"ventral_photophores",
"internal_teeth",
"eye_position",
"oral_gape_shape"
)
) %>%
mutate(trait= gsub("_"," ", trait)) %>%
select(c(trait, depth_layer, mean)) %>%
group_by(
depth_layer, trait
) %>%
summarise(mean_value= mean(mean)) %>%
distinct() %>%
ungroup() %>%
tidyr::pivot_wider(id_cols=depth_layer, values_from = "mean_value", names_from = "trait") %>%
tibble::column_to_rownames(var = "depth_layer")
res.pca <- FactoMineR::PCA(sum_boot_moment, graph = FALSE)
factoextra::fviz_pca_biplot(res.pca, repel = TRUE,
col.var = "gray50",
col.ind = "black",
pointsize = 2,
arrowsize = 0.2,
title = "",
label = "var")+
theme_light()+
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank())
Code
ggsave("PCA_CWM.png", path = "figures/additional_analyses", dpi = 700, height = 5, width = 7)Summarizes bootstrapping output
Extracting raw distributions
Code
# run nonparametric bootstrapping
raw_dist_np <- trait_np_bootstrap(
filled_traits = trait_filling,
raw = TRUE
)
raw_dist_np$depth_layer <- factor(
raw_dist_np$depth_layer,
levels = c(
"Epipelagic",
"Upper mesopelagic",
"Lower mesopelagic",
"Bathypelagic"
)
)
raw_dist_np <- raw_dist_np%>%
mutate(trait= gsub("_"," ", trait)) %>%
filter(!trait%in% c("chin barbel", "fang teeth",
"gland head", "internal teeth",
"large teeth", "retractable teeth",
"small teeth", "ventral photophores"))
ggplot(raw_dist_np, aes(x = log(values), fill = depth_layer, col= depth_layer)) +
geom_density(alpha = 0.1) +
scale_color_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A"))+
scale_fill_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A"))+
labs(x = " log(Trait value)") +
guides(col="none", fill="none")+
facet_wrap(facets = vars(trait), scales = "free")+
theme_light()+
theme(strip.text.x = element_text(size = 10, face = "bold", color="black"),
strip.background = element_rect(fill = "white"),
axis.title = element_text(size = 10),
axis.text = element_text(size = 10))
Functional redundancy
Species Richness (N)
Simpson diversity (D) : probability that two individuals randomly selected from a sample will belong to the same species
Rao diversity (Q) : Sum of distances between pairs of randomly chosen species in trait space weighted by a relative abundance
Functionnal redundancy Ricotta et al. (2016) = 1-Q/D
Code
#Depth biomass data
depth_fish_biomass_indices <- rbind(data_biomass_2002_2019, data_biomass_2021, data_biomass_2022)%>%
as.data.frame()%>%
rename("species"="Nom_Scientifique",
"station"="Code_Station") %>%
left_join(metadata) %>%
select(species, Tot_V_HV, depth, volume_filtered, depth)%>%
# divise biomass by the volume filtered at each trawl (g.m3)
mutate(biomass_cpu=(Tot_V_HV/volume_filtered)*1000)%>%
select(species, depth, biomass_cpu)%>%
replace(is.na(.), 0)%>%
group_by(species, depth)%>%
mutate(biomass=sum(biomass_cpu))%>%
select(-c(biomass_cpu))%>%
distinct()%>%
arrange(depth) %>%
tidyr::pivot_wider(names_from = species, values_from = biomass)%>%
replace(is.na(.), 0)%>%
tibble::column_to_rownames(var = "depth")%>%
#change species name
rename("Lampanyctus_ater"="Nannobrachium_atrum")%>%
rename("Cyclothone_sp"="Cyclothone")%>%
rename("Stomias_boa"="Stomias_boa_boa") %>%
as.matrix()
# Species richness ----
nsp <- as.data.frame(vegan::specnumber(depth_fish_biomass_indices)) %>%
tibble::rownames_to_column(var = "depth") %>%
mutate(depth=as.numeric(depth))
colnames(nsp) <- c("depth", "species_richness")
# Diversity indices
diversity_indices_list <- SYNCSA::rao.diversity(depth_fish_biomass_indices, traits = fish_traits)
diversity_indices <- data.frame(diversity_indices_list[c(2:3)]) %>%
tibble::rownames_to_column(var = "depth") %>%
mutate(depth=as.numeric(depth)) %>%
inner_join(nsp) %>%
mutate(functionnal_redundancy_ricotta= 1-(FunRao/Simpson)) %>%
tidyr::pivot_longer(! depth, names_to = "indices") %>%
# add column with depth layer
mutate(
depth_layer = case_when(
between(depth, 0, 174) ~ "Epipelagic",
between(depth, 175, 699) ~ "Upper mesopelagic",
between(depth, 700, 999) ~ "Lower mesopelagic",
between(depth, 1000, 2000) ~ "Bathypelagic"))
diversity_indices$depth_layer <-
factor(
diversity_indices$depth_layer,
levels = c(
"Epipelagic",
"Upper mesopelagic",
"Lower mesopelagic",
"Bathypelagic"
)
)
diversity_indices$indices <- factor(
diversity_indices$indices,
levels = c(
"species_richness",
"Simpson",
"FunRao",
"functionnal_redundancy_ricotta"
),
labels = c(
"Species richness (N)",
"Simpson diversity (D)",
"Rao diversity (Q)",
"Functional Redundancy R=1-(Q/D)"
)
)
ggplot(diversity_indices, aes(x = depth, y = value)) +
geom_point(alpha = 0.5, size = 1.5) +
geom_smooth(method = "lm", col = "#008E9B", alpha = 0.1) +
facet_wrap(~indices, scales = "free") +
ggpmisc::stat_poly_eq(formula = y ~ x,
aes(label = paste(..rr.label.., ..p.value.label.., ..n.label.., sep = "*`,`~")),
parse = TRUE,
hjust = 0, vjust = 1) +
theme_light() +
theme(axis.title.y.left = element_text(size = 14),
strip.text = element_text(size = 14, face = "bold"),
legend.title = element_text(size = 11),
legend.text = element_text(size = 11),
axis.title.y = element_text(size = 11),
axis.text.y = element_text(size = 12),
strip.background = element_rect(fill = "white"),
strip.text.x = element_text(size = 11, face = "bold", color = "black"))
Code
ggplot(diversity_indices, aes(x =depth_layer, y=value))+
geom_boxplot(aes(col=depth_layer, fill=depth_layer), alpha=0.1) +
scale_fill_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A"))+
scale_color_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A"))+
facet_wrap(~indices, scales = "free")+
theme_light()+
guides(col="none", fill="none")+
theme(axis.text.x = element_blank(),
axis.title.x = element_blank(),
axis.title.y.left = element_text(size =14),
strip.text = element_text(size =14, face="bold"),
legend.title = element_text(size =11),
legend.text = element_text(size =11),
axis.title.y = element_text(size=11),
axis.text.y = element_text(size=12),
strip.background=element_rect(fill="white"),
strip.text.x = element_text(size = 10, face = "bold", color = "black"))
Code
# # Data Preparation
# station_morpho <- morpho_data %>%
# select(station) %>%
# filter(!station%in%c("F0315", "F0367","X0568")) %>%
# arrange(station) %>%
# distinct()
#
# station_sp <- rbind(data_biomass_2002_2019, data_biomass_2021, data_biomass_2022) %>%
# as.data.frame() %>%
# rename("species"="Nom_Scientifique",
# "station"="Code_Station") %>%
# left_join(metadata) %>%
# select(species, Tot_V_HV, volume_filtered, station) %>%
# # divise biomass by the volume filtered at each trawl (g.m3)
# mutate(biomass_cpu=(Tot_V_HV/volume_filtered)*1000) %>%
# select(species, biomass_cpu, station) %>%
# replace(is.na(.), 0) %>%
# group_by(species, station) %>%
# mutate(biomass=sum(biomass_cpu)) %>%
# select(-c(biomass_cpu)) %>%
# distinct() %>%
# tidyr::pivot_wider(names_from = species, values_from = biomass) %>%
# replace(is.na(.), 0) %>%
# arrange(station) %>%
# filter(station %in% station_morpho$station) %>%
# tibble::column_to_rownames(var = "station") %>%
# #change species name
# rename("Lampanyctus_ater"="Nannobrachium_atrum") %>%
# rename("Cyclothone_sp"="Cyclothone") %>%
# rename("Stomias_boa"="Stomias_boa_boa") %>%
# as.matrix()
#
# # Null Model Analysis
# nb_rep <- 1000
#
# resultsRandom_FunRao <-
# matrix(
# NA,
# nrow = nrow(station_sp),
# ncol = nb_rep,
# dimnames = list(rownames(station_sp),
# paste0("Sim.", 1:nb_rep))
# )
#
# resultsRandom_FunRedundancy <-
# matrix(
# NA,
# nrow = nrow(station_sp),
# ncol = nb_rep,
# dimnames = list(rownames(station_sp),
# paste0("Sim.", 1:nb_rep))
# )
#
# for (rep in 1:nb_rep) {
# randomize_mx <- picante::randomizeMatrix(samp = station_sp,
# null.model = "frequency",
# iterations = 100)
#
# sim_cal <- SYNCSA::rao.diversity(randomize_mx , traits = fish_traits)
# sim_FunRao <- sim_cal$FunRao
# sim_FunRedundancy <- sim_cal$FunRedundancy
#
# resultsRandom_FunRao[, rep] <- sim_FunRao
# resultsRandom_FunRedundancy[, rep] <- sim_FunRedundancy
# }
#
# # Comparison with Observed Values
# obs <- SYNCSA::rao.diversity(station_sp, fish_traits)
# obs_FunRao <- obs$FunRao
# obs_FunRedundancy <- obs$FunRedundancy
#
# # Calculate mean and standard deviation of null model FD values for FunRao
# meanNull_FunRao <- rowMeans(resultsRandom_FunRao, na.rm = TRUE)
# sdNull_FunRao <- apply(resultsRandom_FunRao, 1, sd, na.rm = TRUE)
#
# # Calculate mean and standard deviation of null model FD values for FunRedundancy
# meanNull_FunRedundancy <- rowMeans(resultsRandom_FunRedundancy, na.rm = TRUE)
# sdNull_FunRedundancy <- apply(resultsRandom_FunRedundancy, 1, sd, na.rm = TRUE)
#
# # Calculate effect size and standardized effect size for FunRao
# ES_FunRao <- obs_FunRao - meanNull_FunRao
# SES_FunRao <- ES_FunRao / sdNull_FunRao
#
# # Calculate effect size and standardized effect size for FunRedundancy
# ES_FunRedundancy <- obs_FunRedundancy - meanNull_FunRedundancy
# SES_FunRedundancy <- ES_FunRedundancy / sdNull_FunRedundancy
#
# # Combine the results into a data frame
# result_fd_df <- data.frame(
# station = rownames(station_sp),
# SES_FunRao = SES_FunRao,
# SES_FunRedundancy = SES_FunRedundancy
# ) %>%
# inner_join(metadata) %>%
# select(
# -c(
# latitude_start,
# longitude_start,
# latitude_end,
# longitude_end,
# duration,
# distance,
# volume_filtered,
# years,
# comments,
# station
# )
# ) %>%
# mutate(
# depth_layer = case_when(
# between(depth, 0, 174) ~ "Epipelagic",
# between(depth, 175, 699) ~ "Upper mesopelagic",
# between(depth, 700, 999) ~ "Lower mesopelagic",
# between(depth, 1000, 2000) ~ "Bathypelagic"
# )
# ) %>%
# tidyr::pivot_longer(!c(depth, depth_layer), values_to = "values", names_to = "indice")
#
# result_fd_df$depth_layer <- factor(result_fd_df$depth_layer,
# levels = c("Epipelagic", "Upper mesopelagic",
# "Lower mesopelagic", "Bathypelagic"))
#
# # Plot
# ggplot(result_fd_df, aes(x = depth_layer)) +
# geom_point(aes(y = values, col=depth_layer), alpha=0.4) +
# facet_wrap(~indice)+
# geom_boxplot(aes(y = values, col=depth_layer, fill=depth_layer), alpha=0.07) +
# scale_color_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A")) +
# scale_fill_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A")) +
# labs(x="") +
# theme_light()+
# guides(col="none", fill="none")+
# theme(axis.text.x = element_blank(),
# axis.title.x = element_blank(),
# axis.title.y.left = element_text(size =14),
# strip.text = element_text(size =14, face="bold"),
# legend.title = element_text(size =11),
# legend.text = element_text(size =11),
# axis.title.y = element_text(size=11),
# axis.text.y = element_text(size=12),
# strip.background=element_rect(fill="white"),
# strip.text.x = element_text(size = 10, face = "bold", color = "black"))+
# geom_hline(yintercept=0, linetype="dashed", color = "black", linewidth=0.8)SES Functional diversity indices
Standard Effect Size (SES): to eliminate the influence of species richness on the functional diversity indices (Mouchet et al., 2010). Measures the deviation from the random expectation in standard deviation units
null model frequency: Randomize community data matrix abundances (here biomasss) within species (maintains species occurence frequency)
FRic Functional Richness: the proportion of functional space filled by species of the studied assemblage, i.e. the volume inside the convex-hull shaping species. To compute FRic the number of species must be at least higher than the number of functional axis + 1.
FDis Functional Dispersion: the biomass weighted deviation of species traits values from the center of the functional space filled by the assemblage i.e. the biomass-weighted mean distance to the biomass-weighted mean trait values of the assemblage.
FDiv Functional Divergence: the proportion of the biomass supported by the species with the most extreme functional traits i.e. the ones located close to the edge of the convex-hull filled by the assemblage.
FEve Functional Evenness: the regularity of biomass distribution in the functional space using the Minimum Spanning Tree linking all species present in the assemblage.
FOri Functional Originality: the weighted mean distance to the nearest species from the global species pool.
FSpe Functional Specialization: the biomass weighted mean distance to the mean position of species from the global pool (present in all assemblages).
FNND Functional Mean Nearest Neighbour Distance: the weighted distance to the nearest neighbor within the assemblage
FMPD Functional Mean Pairwise Distance: the mean weighted distance between all species pairs.
Code
# Data Preparation ----
station_morpho <- morpho_data %>%
select(station) %>%
filter(!station %in% c("F0315", "F0367", "X0568")) %>%
arrange(station) %>%
distinct()
station_sp <- rbind(data_biomass_2002_2019, data_biomass_2021, data_biomass_2022) %>%
as.data.frame() %>%
rename("species" = "Nom_Scientifique",
"station" = "Code_Station") %>%
left_join(metadata) %>%
select(species, Tot_V_HV, volume_filtered, station) %>%
# Divide biomass by the volume filtered at each trawl (g.m3)
mutate(biomass_cpu = (Tot_V_HV / volume_filtered) * 1000) %>%
select(species, biomass_cpu, station) %>%
replace(is.na(.), 0) %>%
group_by(species, station) %>%
mutate(biomass = sum(biomass_cpu)) %>%
select(-biomass_cpu) %>%
distinct() %>%
tidyr::pivot_wider(names_from = species, values_from = biomass) %>%
replace(is.na(.), 0) %>%
arrange(station) %>%
filter(!station%in%c("H0411","L0731", "L0736")) %>%
tibble::column_to_rownames(var = "station") %>%
# Change species names
rename("Lampanyctus_ater" = "Nannobrachium_atrum") %>%
rename("Cyclothone_sp" = "Cyclothone") %>%
rename("Stomias_boa" = "Stomias_boa_boa") %>%
select(order(colnames(.))) %>%
as.matrix()
# if we want only presence-absence
#station_sp <- replace(station_sp, station_sp > 0, 1)
# dist matrix (gower)
sp_dist_fish <- mFD::funct.dist(
sp_tr = fish_traits,
tr_cat = fish_traits_cat,
metric = "gower",
scale_euclid = "scale_center",
ordinal_var = "classic",
weight_type = "equal",
stop_if_NA = TRUE)
# Compute functional spaces
fspaces_quality_fish <- mFD::quality.fspaces(
sp_dist = sp_dist_fish,
maxdim_pcoa = 10,
deviation_weighting = "absolute",
fdist_scaling = FALSE,
fdendro = "average")
# PC coord
sp_faxes_coord_fish <- fspaces_quality_fish$"details_fspaces"$"sp_pc_coord"
# Calculate functional diversity for the observed data ----
obsFD <- mFD::alpha.fd.multidim(
sp_faxes_coord = sp_faxes_coord_fish[, c("PC1", "PC2", "PC3", "PC4")],
asb_sp_w = station_sp,
scaling = TRUE,
check_input = TRUE,
details_returned = F
)
obsFD_div <- obsFD$functional_diversity_indices
# Null model ----
# Define the number of replications
nb_rep <- 10
# Initialize a list to store results of random functional diversity calculations for each index
indices_names <- colnames(obsFD_div)
resultsRandomFD <- list()
for (index_name in indices_names) {
resultsRandomFD[[index_name]] <- matrix(
NA,
nrow = nrow(station_sp),
ncol = nb_rep,
dimnames = list(rownames(station_sp), paste0("Sim.", 1:nb_rep))
)
}
# Perform randomization and calculate functional diversity for each replication
for (rep in 1:nb_rep) {
randomize_mx <- picante::randomizeMatrix(
samp = station_sp,
null.model = "frequency", # richness, independentswap, trialswap
iterations = 10
)
simFD_cal <- mFD::alpha.fd.multidim(
sp_faxes_coord = sp_faxes_coord_fish[, c("PC1", "PC2", "PC3", "PC4")],
asb_sp_w = randomize_mx,
scaling = TRUE,
check_input = TRUE,
details_returned = F
)
simFD_div <- simFD_cal$functional_diversity_indices
for (index_name in indices_names) {
simFD_index <- simFD_div[, index_name]
# Ensure that simFD_index has the same length as the number of rows in station_sp
if(length(simFD_index) == nrow(station_sp)) {
resultsRandomFD[[index_name]][, rep] <- simFD_index
} else {
stop(paste("The length of", index_name, "does not match the number of rows in station_sp"))
}
}
}
# Initialize dataframes to store mean, standard deviation, effect size, and standardized effect size
meanNullFD <- data.frame(matrix(NA, nrow = nrow(station_sp), ncol = length(indices_names)))
sdNullFD <- data.frame(matrix(NA, nrow = nrow(station_sp), ncol = length(indices_names)))
ES_FD <- data.frame(matrix(NA, nrow = nrow(station_sp), ncol = length(indices_names)))
SES_FD <- data.frame(matrix(NA, nrow = nrow(station_sp), ncol = length(indices_names)))
# Set column names for the dataframes
colnames(meanNullFD) <- indices_names
colnames(sdNullFD) <- indices_names
colnames(ES_FD) <- indices_names
colnames(SES_FD) <- indices_names
# Calculate statistics for each index
for (index_name in indices_names) {
# Calculate mean and standard deviation of null model FD values for each index
meanNullFD[, index_name] <- rowMeans(resultsRandomFD[[index_name]], na.rm = TRUE)
sdNullFD[, index_name] <- apply(resultsRandomFD[[index_name]], 1, sd, na.rm = TRUE)
# Calculate effect size and standardized effect size for each index
ES_FD[, index_name] <- obsFD_div[, index_name] - meanNullFD[, index_name]
SES_FD[, index_name] <- ES_FD[, index_name] / sdNullFD[, index_name]
}
# Combine all results into a single dataframe
results_df <- cbind(
# obsFD_div,
# meanNullFD = meanNullFD,
# sdNullFD = sdNullFD,
# ES_FD = ES_FD,
SES_FD = SES_FD
)
# Add row names
rownames(results_df) <- rownames(station_sp)
# Plot ----
# Output the results dataframe
results_df_plot <- results_df %>%
tibble::rownames_to_column(var = "station") %>%
inner_join(metadata %>% select(station, depth), by = "station") %>%
mutate(
depth_layer = case_when(
between(depth, 0, 174) ~ "Epipelagic",
between(depth, 175, 699) ~ "Upper mesopelagic",
between(depth, 700, 999) ~ "Lower mesopelagic",
between(depth, 1000, 2000) ~ "Bathypelagic"
)
) %>%
select(-station) %>%
tidyr::pivot_longer(!c(depth, depth_layer),
values_to = "values",
names_to = "indice") %>%
mutate(indice = stringr::str_replace(indice, "^SES_", "")) %>%
filter(
indice %in% c(
"FD.fric",
"FD.fdis",
"FD.fmpd",
"FD.fnnd",
"FD.feve",
"FD.fdiv",
"FD.fori",
"FD.fspe"
)
)
results_df_plot$depth_layer <- factor(results_df_plot$depth_layer,
levels = c("Epipelagic", "Upper mesopelagic",
"Lower mesopelagic", "Bathypelagic"))
results_df_plot$indice <- factor(results_df_plot$indice,
levels = c("FD.fric",
"FD.fdis",
"FD.fdiv",
"FD.feve",
"FD.fori",
"FD.fspe",
"FD.fnnd",
"FD.fmpd"
))
ggplot(results_df_plot, aes(x = depth_layer, y = values, fill = depth_layer)) +
facet_wrap(~indice, scales = "free") +
geom_point(alpha = 0.5, size = 1, position = position_jitter(width = 0.2), aes(col=depth_layer)) +
geom_boxplot(alpha = 0.1, outlier.shape = NA, width = 0.5, aes(col=depth_layer, fill=depth_layer))+
scale_color_manual(values = c("#FEA520", "#D62246", "#6255B4", "#3C685A")) +
scale_fill_manual(values = c("#FEA520", "#D62246", "#6255B4", "#3C685A")) +
labs(
x = "",
y = "Standard Effect Size (SES)"
) +
theme_light() +
theme(axis.text.x = element_blank(),
strip.text.x = element_text(
size = 13,
face = "bold",
color = "gray20"
),
strip.background = element_rect(fill = "white"),
axis.title = element_text(size = 13),
axis.text = element_text(size = 13))+
geom_hline(yintercept = 0, linetype = "dashed", color = "black", linewidth = 0.8)+
guides(col="none", fill="none")
Isotopes data
- 4 sp manquantes en d15N = P. coregonoides, E. pelacanoides, H. anomala, S. bathyphilus
- functional space : colours according to d15N values
- significant relationships between d15n and traits
- categorical traits
Plan
I. Key traits and foraging strategies by depth layer
Figure 1: Lateral-view photographs
Figure 2: Functional spaces
At the community level & by depth layer
Only PC1-PC2 and PC3-PC4 in appendices?
Table 1: significant traits
which traits significantly influence each axis?
at the community level & by depth layer
Figure 3: Community-Weighted Mean
II. Functional diversity indices
Figure 4: Standardized Effect Size
- Trait convergence in the epipelagic layer and trait divergence in the bathypelagic layer
III. Functional rarity
Figure 5: Uniqueness-geographical restrectiveness
